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The  finite-element  technique  has  the  potential  to  provide  a  very  accurate  treatment  of  the 
physics  of  acoustic-wave  propagation  in  inhomogeneous  media.  This  article  describes  the 
development  of  a  finite-element  model  for  acoustic  propagation  in  complex  ocean 
environments  and  its  validation.  The  computational  model  can  handle  range  and  depth 
dependence  in  both  sound  speed  and  density,  as  well  as  rapid  variations  in  bottom  topography. 
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INTRODUCTION 

There  are  many  underwater  acoustic  propagation  mod¬ 
els  currently  in  use  in  the  scientific  community,  each  of 
which  has  its  own  particular  advantages  and  disadvantages. 
The  need  for  more  accurate  models,  particularly  for  applica¬ 
tions  at  low  frequencies  where  there  can  be  significant 
boundary  penetration,  is  growing.  It  is  essential  to  give  a 
more  accurate  treatment  of  the  physics  of  acoustic/elastic- 
wave  propagation.  This,  in  turn,  means  that  the  computa¬ 
tional  models  must  include  the  effects  of  depth-  and  range- 
dependent  sound  speeds,  bottom  topography,  and 
range-dependent  shear  in  ocean  sediments.  The  finite-ele¬ 
ment  method  has  the  potential  to  provide  a  very  accurate 
treatment  of  the  physics  of  wave  propagation  in  such  com¬ 
plex  ocean  environments. 

The  use  of  finite-element  techniques  to  model  acoustic 
waves  is  not  new.  Kalinowski  has  provided  an  excellent  sum¬ 
mary  of  the  work  of  many  authors,  including  his  own,  on  the 
application  of  finite  elements  to  various  acoustic  problems. 1 
Earlier  work  has  demonstrated  that  finite  elements  can  be 
used  to  describe  the  effect  of  bottom  shear  on  propagation 
over  distances  of  a  few  wavelengths.  A  judicious  blending  of 
finite  elements  and  the  boundary  integral  method  has  been 
used  to  study  sound-structure  interactions  for  submerged 
bodies.  Kuo  and  co-workers  have  used  the  finite-element 
method  to  model  the  time  domain  pulse  propagation  ob¬ 
served  in  laboratory  scale  models.2  The  two-dimensional  fi¬ 
nite-difference  model  of  Stephen  is  similar  in  that  it  models 
pulse  propagation  and  includes  the  effect  of  shear  in  the 
ocean  bottom. 1  These  models  are  well  suited  to  study  various 
types  of  geophysical  phenomena  in  the  time  domain.  Gold¬ 
stein  etal.  have  studied  accelerated  convergence  methods  for 
iterative  numerical  solution  techniques  applicable  to  the 
large  systems  of  equations  associated  with  finite-element 
models.4  These  iterative  techniques  are  especially  important 
when  one  deals  with  numerical  solutions  of  elliptic  equations 
that  require  large  quantities  of  computer  memory. 

This  article  describes  a  finite-element  ocean  acoustic 
model  (FOAM)  which,  in  its  present  form,  can  very  accu¬ 


rately  model  propagation  over  several  kilometers.  It  handles 
depth-  and  range-dependent  sound  speed,  density,  and  al¬ 
most  arbitrary  variations  in  bottom  topography.  It  does  not 
yet  include  shear  effects,  and  compared  to  some  other  types 
of  models  such  as  the  parabolic  equation  (PE)  model,  it  is 
computationally  slow.  However,  unlike  the  one-way  PE 
models  and  most  normal-mode  models,  it  is  a  full-wave  solu¬ 
tion,  i.e.,  it  includes  backscatter  and  both  the  discrete  modes 
and  continuous  spectra.  There  are  no  restrictions  on  angles 
of  propagation.  Unlike  the  time-domain  models  previously 
mentioned,  which  are  often  applied  to  seismic  phenomena, 
ours  is  a  frequency-domain  model  designed  primarily  for 
transmission  loss  studies. 

In  Sec.  I,  the  finite-element  mesh  utilized  in  the  model  is 
discussed,  as  is  the  particular  form  of  the  basis  or  interpola¬ 
tion  functions.  Then,  in  Sec.  II,  the  development  of  the  finite- 
element  equations  is  sketched,  followed  by  the  numerical 
solution  technique  in  Sec.  III.  Two  of  the  ways  in  which  the 
physics  and  the  computer  code  were  verified  are  described  in 
Sec.  IV.  A  brief  summary  of  our  results  and  conclusions,  as 
well  as  future  plans  to  continue  the  development  of  the  mod¬ 
el  are  given  in  Sec.  V. 


I.  FINITE-ELEMENT  MESH  AND  INTERPOLATION 
FUNCTIONS 

In  the  application  of  the  finite-element  technique  to 
solving  partial  differential  equations,  the  differential  equa¬ 
tion,  and  its  unknown  solution  are  replaced  by  a  system  of 
algebraic  equations  in  terms  of  the  parameters  defining  an 
approximate  solution.  One  partitions  the  domain  of  the 
problem  into  nonoverlapping  elements  and  assumes  a  simple 
form  for  the  approximate  solution  within  each  element. 
These  local  representations  are  then  joined  together  by  using 
the  appropriate  physical  continuity  conditions  to  provide  a 
global  solution.  The  partial  differential  equations  and  asso¬ 
ciated  boundary  conditions  that  arise  in  many  science  and 
engineering  problems  give  matrix  systems  of  equations 
which  are  sparse,  banded,  and  often  symmetric.  Fortunate- 


1478 


J.  Acoust  Soc  Am  86  (4),  October  1989 


0001 -4966/89/1 01 478-06$00  80 


©  1 989  Acoustical  Society  of  America 


1478 


ly,  that  is  the  case,  too,  for  the  description  of  ocean  acoustic- 
wave  propagation.  The  reader  unfamiliar  with  finite-ele¬ 
ment  techniques  may  want  to  consult  the  excellent  text  by 
Reddy.5 

Consider  the  typical  problem  of  ocean  acoustics  sug¬ 
gested  by  Fig.  1.  An  acoustic  source  is  placed  somewhere  in 
the  water  column.  The  sound  speed,  whether  in  the  water  or 
in  the  ocean  bottom,  may  vary  with  both  depth  and  range. 
The  depth  of  the  water-bottom  interface  may  vary  with 
range.  The  acoustic  pressure  is  taken  to  vanish  at  the  air- 
water  interface.  If  the  ocean  bottom  is  assumed  to  be  homo¬ 
geneous  beyond  some  maximum  depth,  and  if  the  sound- 
speed  profiles  do  not  vary  beyond  some  maximum  range, 
then  along  the  right  and  bottom  sides  of  the  domain  to  be 
modeled  one  needs  to  impose  some  kind  of  approximate 
boundary  condition  which  does  not  reflect  any  energy  back 
into  the  interior. 

Although  basis  or  interpolation  functions  can  be  con¬ 
structed  for  many  types  of  polygonal  regions,  the  triangle  is 
one  of  the  most  widely  used  finite  elements.  At  present,  the 
computer  model  FOAM  makes  use  of  linear  Lagrange  inter¬ 
polation  functions.  For  each  finite  element  ft‘‘  that  makes  up 
the  domain  to  be  modeled,  one  must  form  a  set  of  interpola¬ 
tion  funcf'^ns.  The  unknown  solution  will  be  expressed  as  a 
linear  c  ination  of  these  interpolation  functions.  In  addi¬ 
tion  to  vanishing  outside  the  element  Slr,  these  functions  sat¬ 
isfy  the  following  conditions: 


< ifi'/'UjS,)  -  8jj, 

(1) 

t I’l'Hz.r)  =  1, 

(2) 

<  -  i 


where  r,,z,  is  the  location  of  the  /th  vertex  of  the  element. 
Finite-element  models  generally  make  use  of  two  node  num- 
►  bering  schemes:  a  local  one  as  used  here,  i—  1,  2,  or  3,  and  a 
global  one,  I  =  1,2,...,«,  where  n  is  the  total  number  of  nodes 
in  the  mesh.  Given  the  element  number,  one  can  map  from 
the  local  node  numbers  to  the  global  node  numbers. 

For  triangular  elements  with  linear  interpolation  func¬ 
tions,  the  «/’,  are  given  by 

tl'Y'(z.r)  =  (a,  +0,z  +  r,r)/2Ac,  (3) 

where  a,,  0,,  and  y,  are  the  constants 

a,=zlrk-zkrr  0,  =  r, -rA,  y,=zk-zr  (4) 


FIG.  1.  Schematic  drawing  illustrating  the  finite-element  mesh,  boundary 
conditions,  and  nodes  adjusted  to  follow  the  ocean  bottom  bathymetry. 


In  the  above  equations,  i^j^k,  and  /,/.  and  k  permute  in  a 
natural  order.  The  constant  At,  is  the  area  of  the  triangle. 

Integrals  involving  products  of  two  or  more  interpola¬ 
tion  functions  can  be  simplified  by  the  use  of  area  coordi¬ 
nates  Z.,,  where 

L,(z,r)  =  A,/Ae  =  tp,(z,r).  (5) 

Geometrically,  A,  is  the  area  of  the  triangle  formed  by  the 
point  ( z,r )  and  nodesy  and  k.  The  following  integration  for¬ 
mula  is  quite  useful  for  evaluating  several  of  the  required 
integrals: 

f  (£,)'"(/., )"(£., VdA  = - - 2 A,  ,  (6) 

J  (m  +  n  +  p  +  2 ) ! 

where  m,  n,  and  p  are  arbitrary  non-negative  integers. 

In  our  computer  model,  we  have  used  triangular  three- 
node  elements,  as  previously  discussed,  so  that  the  pressure 
within  each  element  will  be  expressed  as  a  linear  combina¬ 
tion  of  these  interpolation  functions: 

3 

Pu’(z,r)  =  £  PjW'izs).  (7) 

;  -  i 

Extensions  to  four-node,  six-node,  etc.  triangles  will  allow 
more  accurate  representations  of  the  pressure  field  since 
these  require  quadratic,  cubic,  and  higher-order  polynomial 
expressions  for  the  interpolation  functions. 


II.  THEORY 

The  time-harmonic  complex  pressure  resulting  from  a 
point  source  of  unit  strength  and  angular  frequency  w  at 
position  r,  in  an  inhomogeneous  ocean  is  determined  by  the 
wave  equation 

pV- (p-'VP)  +k2P=  -<5(r-rJ,  (8) 

and  appropriate  boundary  conditions.  We  will  make  the 
common  assumption  of  an  ocean  environment  axially  sym¬ 
metric  about  a  vertical  line  passing  through  the  source,  so 
that  in  the  above  equation,  p  =  p(z,r)  is  the  density, 
k  —  k(z,r)  =  co/c(z,r)  is  the  wavenumber  and  cU,r )  is  the 
sound  speed  To  allow  for  attenuation,  the  wavenumber  can 
be  given  an  imaginary  part, 

k  =  (co/ c) ( 1  -t-  irja),  (9) 

where  a  is  the  attenuation  of  the  medium  in  units  of  dB  per 
wavelength  and  rj  =  {40n  log,,,  e)  - '. 

The  finite-element  technique  will  eventually  yield  a  set 
of  pressure  values  at  the  nodes  of  the  grid  as  illustrated  in 
Fig.  1 .  One  of  the  key  features  of  the  method  is  the  ability  to 
use  a  grid  of  elements  adapted  to  the  specific  problem.  For 
example,  if  the  density  has  a  discontinuous  jump  in  going 
across  the  water-bottom  interface,  the  mesh  would  be 
formed  so  that  the  interface  coincides  with  the  boundaries  of 
the  appropriate  interior  elements. 

A.  Finite-element  equations 

One  uses  the  wave  equation  within  each  element  (and 
the  interpolation  functions  defined  for  that  element)  to  ar¬ 
rive  at  its  contribution  to  the  system  of  equations.  Let  us 
consider  an  element  IT  that  does  not  contain  the  source. 
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Multiply  Eq.  (8)  by  a  test  function  V and  integrate  over  the 
element  O'  : 


VPrdA 


' k'-vPrdA 


where  the  divergence  theorem  has  been  used  to  transform 
one  of  the  integrals.  The  integrals  are  over  the  volume  of  the 
element  iV  and  around  the  boundary  P  of  the  element.  Be¬ 
cause  of  the  assumption  of  axisymmetry,  the  solution  and 
the  test  function  are  taken  to  depend  upon  z  and  r  only.  The 
volume  differential  is  dr  =  2irr  dr  dz  —  2rr  dA,  and  a  com¬ 
mon  factor  of  2r  has  been  canceled  everywhere.  The  partial 
derivative  in  the  last  in‘egrand  is  with  respect  to  the  normal 
to  the  element's  boundary. 

Take  the  test  function  v  to  be  equal  to  one  of  the  interpo¬ 
lation  functions,  say,  tA, ,  and  in  the  first  two  integrals  express 
the  pressure  within  the  element  in  the  form 

p=  £  p, (id 


Eq.  (10)  becomes, 


i(J 

i  t  \Jir 


Vil^rdA 


'k  H’Y  'il'l'-'rdA  ^Pj 


If  the  density  p  and  wavenumber  k  are  treated  as  constants 
within  each  element,  the  first  and  second  integrals  in  Eq. 
(12)  can  be  evaluated  exactly.  They  are  proportional  to  the 
following  two  integrals: 


7’  =  J 

t  r = | 


Vtty’rdA, 


The  second  of  these  can  be  evaluated  using  the  substitution 
r=  X  rkA  (15) 

k  I 

and  then  the  integration  formula  given  above  in  Eq.  (6). 
Note  that  both  K  b’*  and  M  ]/’  are  symmetric  matrices.  They 
are  usually  referred  to  as  stiffness  and  mass  matrices,  respec¬ 
tively. 


B.  Boundary  integrals  and  boundary  conditions 

For  those  elements  that  are  in  the  interior,  the  boundary 
integral  is  not  zero,  but  its  contributions  to  the  global  system 
of  equations  will  exactly  cancel  with  like  terms  coming  from 
neighboring  elements.  One  only  need  recall  that  the  term 
p  'dP  /dn  is  proportional  to  particle  displacement  normal  to 
the  boundary  and  it  must  be  continuous  across  the  bound¬ 
ary.  For  the  case  where  the  element  has  one  or  two  of  its  sides 
on  an  actual  boundary  of  the  FE  model,  one  of  three  possible 
boundary  conditions  are  allowed:  ( 1 )  Neumann  boundary 


condition,  dP  /dn  =0,  (2)  Dirichlet  boundary  condition, 
/*  =  0,  or  (3)  radiating  boundary  condition.  The  first  two 
conditions  are  treated  exactly  in  the  computer  model  using 
standard  finite-element  techniques.  In  the  first  case,  the 
boundary  integral  is  zero:  so  the  first  two  integrals  discussed 
are  the  only  ones  that  contribute  anything.  In  the  second 
case,  the  equations  involving  the  nodes  for  which  P  =  0  are 
simply  replaced  with  P  =  0.  The  third  can  only  be  treated 
approximately  and  requires  considerably  more  work. 

The  radiating  boundary  condition  means  that  waves 
should  pass  through  the  boundary  without  reflection.  Ideal¬ 
ly.  one  imposes  a  radiation  condition  infinitely  far  from  the 
source.  Usually,  in  a  computer  model  one  must  use  an  ap¬ 
proximate  radiation  condition  on  a  finite  boundary.  At  pres¬ 
ent,  in  FOAM,  this  can  be  treated  approximately  by  one  of 
two  boundary  conditions — a  narrow-angle  radiation  bound¬ 
ary  condition  or  a  wide-angle  radiation  boundary  condition. 
These  are  both  described  below.  The  narrow-angle  radiation 
condition,  so-called  because  it  is  valid  for  plane  waves  inci¬ 
dent  on  the  boundary  within  a  small  angular  interval  about 
the  normal,  is  equivalent  to  a  damping  condition  used  earlier 
by  Kalinowski.1  It  is  given  by  the  equation 

dP 

— —  =  ikP.  (16) 

dn 

The  wide-angle  radiation  condition  that  we  have  developed 
is  useful  for  a  somewhat  wider  interval  about  the  normal 
direction  and  is  actually  a  kind  of  PE  used  as  a  boundary 
condition: 


*l  =  ikP+n\d!Pt 
dn  V  2k  )  dr 


where  the  derivatives  d  2/ds 2  are  taken  tangent  to  the  bound¬ 
ary.  A  comparison  that  shows  the  advantage  of  the  second 
radiating  boundary  condition  over  the  first  is  illustrated  in 
Fig.  2(a),  where  the  effective  reflection  coefficients  for  these 
two  boundary  conditions  are  plotted.  The  percent  of  inci¬ 
dent  energy  reflected  is  shown  in  Fig.  2(b).  One  can  con¬ 
clude  that  the  wide-angle  radiation  condition  allows  a  higher 
incident  angle  to  be  reached  before  the  solution  is  seriously 
contaminated  by  false  reflections. 


100  R I 
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FIG  2  Comparisons  of  (a)  reflection  coefficients  and  (h)  intensities  as  a 
tune! loll  of  incident  angle  for  the  narrow-angle  (NAi  and  wide-angle  1®® 
(WA|  radiation  conditions.  ,r 
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Upon  substituting  these  approximate  radiating  bound¬ 
ary  conditions  into  Eq.  (12)  in  the  last  term,  one  finds  that 
their  contributions  can  be  written  in  terms  of  one  or  both  of 
the  integrals: 


Both  are  symmetric  matrices,  and  the  global  equations  re¬ 
sulting  from  either  of  these  boundary  conditions  are  still 
symmetric,  an  important  feature  for  the  solution  technique. 
A  small  amount  of  attenuation  can  be  introduced  near  any 
radiating  boundaries  to  further  extinguish  unwanted  reflec¬ 
tions. 

Once  the  contributions  of  the  various  elements  is  deter¬ 
mined,  the  global  system  of  equations  is  formed  by  mapping 
the  local  node  numbers  onto  the  global  node  numbers,  giving 
rise  to  the  global  pressure  vector  P,  and  combining  all  of  the 
subsystems  into  a  single  global  system  of  equations.  The 
source  contribution  appears  as  a  single  entry  in  the  global 
force  vector  GF,  and  the  stiffness,  mass,  and  boundary  con¬ 
tributions  finally  yield  a  global  stiffness  matrix  GSTIFF. 
One  then  must  solve  the  system 

GSTIFF  •  P  =  GF,  (20) 

where  P  and  GF  are  ^-dimensional  vectors,  and  GSTIFF  is 
the  n  X  n  system  matrix;  here,  n  is  the  total  number  of  pres¬ 
sure  nodes  in  the  finite  element  mesh. 

III.  NUMERICAL  SOLUTION 

The  global  stiffness  matrix  is  a  complex,  symmetric  ma¬ 
trix  that  can  be  factored  in  the  following  form'’: 

GSTIFF  =  UtDU,  (21) 

where  U  is  an  upper  triangular  nXn  matrix  with  ones  along 
its  diagonal,  and  D  is  an  n  X  n  diagonal  matrix.  The  most 
time  consuming  part  of  the  calculations  is  the  factorization 
of  GSTIFF.  One  can  generate  the  £/,,  and  D 'u  in  a  column¬ 
wise  fashion  by  using  the  recurrence  equations  given  below 
for  j  =  2,3 

v„  =  (^-)  (X  -  l’  Dkk  Uk,  £4,),  1  <  /  <j,  ( 22 ) 

DU=An->l  V  Ij,  l  <i  —  j.  (23) 

A  1 

Afterwards,  the  complex  pressure  is  obtained  in  three  steps. 
In  the  first  step,  one  solves  for  the  complex  vector  Z  in  the 
lower  triangular  system, 

UTZ  =  GF.  (24) 

Since  U1  is  a  lower  triangular  matrix,  the  Z,  can  be  calculat¬ 
ed  in  a  series  of  forward  substitutions: 

Z,  =  GF,  -  £  UllGFr  />  1.  (25) 

/  i 

The  second  step  consists  of  solving  for  the  complex  vector  Y 
defined  by  the  diagonal  system 

DY  =  Z.  (26) 


Since  D  is  a  diagonal  matrix,  the  Y,  are  simply  equal  to 
Z,/D, .  Finally,  the  desired  pressure  vector  results  from  solv¬ 
ing  the  upper  triangular  system: 

UP  =  Y.  (27) 

The  elements  of  P  are  found  by  a  series  of  back  substitutions: 

P,  =  Y,  —  £  U,JPJ,  i  <n.  (28) 

J  -  i  *  l 

The  computer  model  uses  a  banded  version  of  a  factori¬ 
zation  algorithm  to  take  full  advantage  of  the  banded  nature 
of  GSTIFF  and  U  to  minimize  the  storage  requirements.7 

IV.  MODEL  VERIFICATION 

In  addition  to  solving  a  number  of  boundary  value  prob¬ 
lems  with  boundary  conditions  selected  to  yield  solutions 
with  particular  types  of  symmetries,  the  model  has  been  used 
to  solve  two  problems  more  familiar  to  the  underwater 
acoustics  community,  i.e.,  the  Lloyd’s  mirror  effect  and  an 
upslope  wedge  propagation  benchmark  problem  recently 
discussed  at  Acoustical  Society  of  America  meetings. 

A.  Lloyd’s  mirror  effect 

As  an  early  test  of  the  accuracy  of  our  FE  model,  the 
Lloyd’s  mirror  effect  was  examined.  This  image-interference 
effect  occurs  in  underwater  acoustics  when  the  acoustic 
source  is  near  the  sea  surface  and  the  sea  surface  is  not  very 
rough  (i.e.,  the  sea  surface  acts  as  a  pressure-release  sur¬ 
face).  Then,  an  interference  pattern  in  the  sound  field  results 
from  constructive  and  destructive  interference  between  the 
direct  and  surface-reflected  acoustic  waves.  There  is  an  ana¬ 
lytic  solution*  that  was  used  as  a  baseline  result  against 
which  the  FOAM  calculations  could  be  compared.  Figure  3 
shows  two  contour  plots  of  constant  pressure  at  2-dB  inter¬ 
vals  for  an  upper  half-space  of  air  and  a  lower  half-space  of 


RANGE  (KM) 


FIG.  .V  Contour  plots  of  the  I, loul's  mirror  effect  computed  b>  ( a )  an  exact 
solution  and  (h)  FOAM  using  the  wide-angle  radiation  condition  Con¬ 
tours  start  at  0  dt)  and  increase  in  intervals  of  d  dB. 
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water  with  0.5  dB  per  wavelength  absorption.  The  25-Hz  cw 
source  was  placed  100  m  below  the  surface.  The  familiar 
Lloyd’s  mirror  beams9  are  quite  evident.  Unlike  the  analytic 
solution,  our  model  must  simulate  the  lower  half-space  using 
a  bounded  medium  terminating  with  a  radiation  boundary 
condition.  The  results  shown  in  Fig.  3(b)  were  obtained  us¬ 
ing  the  wide-angle  radiation  boundary  condition  on  the  bot¬ 
tom  and  right  boundaries  of  our  model,  and  compare  favor¬ 
ably  with  the  analytic  solution  [Fig.  3(a)]  except  near  the 
radiating  boundaries.  Some  low  grazing  angle  reflections  are 
evident  along  the  bottom  and  right  side  as  can  be  seen  by 
comparing  the  two  figures  below. 

Although  this  appears  to  be  a  very  simple  problem,  i.e., 
a  point  source  in  a  homogeneous  half-space  with  a  pressure 
release  surface,  it  is  one  that  cannot  be  done  very  well  with 
most  PE’s  since  they  cannot  handle  the  steep  beams,  nor  will 
they  give  correct  results  at  very  short  ranges.  It  also  provides 
a  rather  stringent  test  of  the  radiation  boundary  condition 
that  we  have  developed  for  the  bottom  and  far  right  boun¬ 
daries  of  the  computational  model. 

B.  ASA  benchmark  problem 

In  1987  at  the  1 13th  meeting  of  the  Acoustical  Society  of 
America  (ASA),  a  special  session  was  devoted  to  numerical 
solutions  of  several  types  of  benchmark  problems.10  One 
problem  dealt  with  upslope  propagation  in  a  wedge-shaped 
underwater  channel.  A  25-Hz  cw  acoustic  source  was  placed 
at  a  depth  of  1  'X)  m  below  the  surface  of  the  water.  The  ocean 
bottom  was  initially  200  m  deep  and  gradually  decreased  in 
depth  at  a  2.86-deg  rise  until  it  intersected  the  ocean  surface 
at  a  range  of  4.0  km.  The  density  of  the  water  was  1  g/cc  and 
the  sound  speed  was  1500  m/s.  The  wedge  half-space  had  a 
density  ofl.5  g/cc,  a  sound  speed  of  1700  m/s,  and  an  atten- 
utiou  of  0.5-dB/wavelength.  In  addition  to  the  continuous 
spectia,  the  source  excites  three  propagating  modes  that  pass 
through  their  cutoff  depths  at  ranges  of  approximately  813, 
2088,  and  3363  m. 

A  range-dependent  coupled-mode  model"  was  used  as 
the  benchmark  result12  against  which  other  models  were 
compared.  Figure  4  shows  a  comparison  of  our  FOAM  re¬ 
sults  (solid  curve)  compared  to  the  coupled-mode  bench¬ 
mark  model  results  (broken  curve).  The  transmission  loss 
curves  are  almost  indistinguishable.  A  contour  plot  showing 
transmission  loss  is  given  in  Fig.  5.  One  can  see  the  effect  of 
mode  dumping  at  the  appropriate  ranges.  Both  of  these  mod¬ 
els  are  full-wave  range-dependent  models  and  can  account 
for  the  backscatter  and  interference  that  comes  from  the 
wedge.  Clearly,  the  finite-element  results  compare  very  fa¬ 
vorably  with  the  benchmark  results.  Because  the  many  pre¬ 
senters  at  the  benchmark  session  used  a  number  of  different 
computers  and  did  not  all  give  epu  times,  we  are  not  able  to 
make  a  comparison  of  relative  computational  effort. 

V.  SUMMARY  AND  CONCLUSIONS 

The  development  of  a  finite-element  model  for  acoustic- 
propagation  in  complex  ocean  environments,  FOAM,  has 
been  described.  Results  of  two  tests  used  to  verify  the  model, 
the  Lloyd’s  mirror  effect,  and  an  upslope  wedge  propagation 
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FIG.  4.  Transmission  loss  plots  (in  dB)  for  the  ASA  benchmark  wedge 
problem,  comparing  FOAM  with  the  ASA  benchmark  solution  obtained 
from  the  Evans'  coupled  mode  model. 


problem  discussed  extensively  at  recent  ASA  meetings  have 
been  shown.  The  transmission  loss  calculations  using  the 
finite-element  model  compare  very  favorably  with  the  analy¬ 
tical  and  benchmark  solutions.  Because  of  its  potential  to 
provide  a  very  complete  treatment  of  propagation  in  com¬ 
plex  ocean  environments,  we  believe  that  it  is  important  to 
continue  its  development,  particularly  to  include  shear  ef¬ 
fects  in  the  ocean  bottom  and  to  extend  it  to  three-dimen¬ 
sional  environments,  while  at  the  same  time  working  to  im¬ 
prove  its  computational  efficiency  and  to  remove  its  current 
restriction  to  short  ranges.  The  latter  is  principally  a  com¬ 
puter  memory  problem,  not  a  restriction  of  the  finite-ele¬ 
ment  technique  itself. 


FIG  5.  Contour  plot  ( in  dB)  from  FOAM  for  the  ASA  benchmark  wedge 
problem 
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